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ABSTRACT 

This paper presents a comparison of the predictions for the 2 and 3-point correlation 
functions of density fluctuations, ^ and C, in gravitational perturbation theory (PT) 
against large Cold Dark Matter (CDM) simulations. This comparison is made possible 
for the first time on large weakly non-linear scales (> 10/i“^Mpc) thanks to the devel¬ 
opment of a new algorithm to estimate correlation functions for millions of points in 
only a few minutes. Previous studies in the literature comparing the PT predictions 
of the 3-point statistics to simulations have focused mostly on Fourier space, angular 
space or smoothed fields. Results in configuration space, such as the ones presented 
here, were limited to small scales were leading order PT gives a poor approximation. 
Here we also propose and apply a method for separating the first and subsequent orders 
contributions to PT by combining different output times from the evolved simulations. 
We find that in all cases there is a regime were simulations do reproduce the leading 
order (tree-level) predictions of PT for the reduced 3-point function ~ • For 

steeply decreasing correlations (such as the standard CDM model) deviations from 
the tree-level results are important even at relatively large scales, cs 20 Mpc/h. On 
larger scales ^ goes to zero and the results are dominated by sampling errors. In more 
realistic models (like the ACDM cosmology) deviations from the leading order PT 
become important at smaller scales r ~ 10 Mpc/h, although this depends on the par¬ 
ticular 3-point configuration. We characterise the range of validity of this agreement 
and show the behaviour of the next order (1-loop) corrections. 


1 INTRODUCTION 

The systematic use of the 2-point and 3-point density cor¬ 
relation functions in the context of the large scale structure 
in the universe was introduced by Pe ebles and his collabo¬ 
rators in the be ginning of the 70’s ([Peebles 1973 ; Hauser 
fc_j)ggbjes_jll7^ Peebles fc Hauser 1974 Peebles 19^ ; Pee- 
bles fc Groth 197f| ). It is well known that a Gaussian field 
produces a zero 3-point correlation function. In the stan¬ 
dard cosmological model, the simplest (1-field) inflationary 
scenarios find a Gaussian (or almost Gaussian) form for 
the distribution of initial matter density and fluctuations 


^ It is therefore important to study if there is the pos¬ 
sibility of disentangling the non-gaussian effects produced 
by gravity from the non-gaussian characteristics of initial 
conditions. We would also like to use the 3-point func¬ 
tion to differentiate among different cosmological models 
through their distinctive (non-linear) gravitational evolu¬ 
tion. As mentioned above, during the last 25 years there have 
been many works related to the study of the 3-point corre¬ 
lation function. Most of the recent studies, have centered 


on the projected distribution (Gaztanaga & Frieman 199£; 


& Kaniiuukuwski 2000 y Gangui Mai Lin 200C ). lluwevei, 
there are more complex inflationary models (Allen et al. 
19^ Kofonan^^PogosyanM98^ salopek et al. 1989; [Linde 


Buchalter fc Kam ionkowski 2000 1) or its Fourier counterpart, 
the bispectrum (pcoccimarro fc Frieman 1996; [Scoccimarro 


I-1 rr -1 IT*-1 I-1 uisptJULi uiii (pcucuiiiicniu 06 nieiiiciii uyyu, pcuuciiiiciiiu 

(iFalk et al. 199^; [GaiigUl et al. 1994^ iGaligUl 1^ [^^/allg| IScoccimarro 20001). In this work we deal lith the study 


& Mukhanov 1997; Peebles_^999ajb and also models based on 
topologi c al defects (T Vllenki^^^: ^ Hill et 

al. 1986; Furok 1989r Albrecht fc Stebbins 1992 Gaztanaga 
fc Mahoenen 1996 ) which predict non-gaussian initial con¬ 
ditions for primordial perturbations. 

Even if we start from gaussian initial conditions, gravity 
makes fluctuations evolv e into a non-gaussian distribution 
(Bernardeau et al. 2001) as will also be shown in section 


of how to compute the 3-point correlation function from sim¬ 
ulations in real 3-D (or configuration) space and its com- 
parison to PT. This topi c has been investigat ed before by 
(Matsubara fc Suto 1994; ling fc Borner 1997) but, due to 


difficulties in dealing with such a large number of particles 
of the simulations, the study was only done with random 
subsamples and only at small scales. It remains to be seen 
if this is the reason why in these previous attempts it was 
found that N-body simulations do not reproduce well the 
leading order PT predictions. 
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2 J. Barriga & E. Gaztanaga 


In this work we develop a means of dealing with such 
large numbers, explained in section and also a method 
to extract non-linearities from simulations in order to com¬ 
pare the result more properly to perturbation theory values, 
in section 3.3. The actual comparison with simulations is 
presented in section §4. 


2 THE 2 AND 3-POINT FUNCTION 

In this section we recall the definitions of the 2 and 3-point 
functions. We also revised the predictions for these quan¬ 
tities according to gravitational evolution in perturbation 
theory starting from gaussian initial conditions. 


2.1 Definitions 


The 2-point correlation function ^{ri 2 ) is defined by the joint 
probability of finding one galaxy (or a point) in a small vol¬ 
ume SVi and another one in the volume SV 2 separated by 
a distance ri 2 , when choosing the two volumes randomly 
within a large (that is, representative) volume of the uni¬ 
verse, 

5P2 = p" [1 + e(ri2)] 5 V 1 SV 2 (1) 

where p is the mean number of galaxies per unit volume. In 
other words, 5 (ri 2 ) represents the excess probability, com¬ 
pared to a random distribution, of finding a galaxy at a 
distance ri 2 from another galaxy. When we approximate 
the distribution of objects by a continuous density function 
p(r) = p[l -I- (5(f^] we can write ^ in terms of density fluctu¬ 
ations: 


^(^■ 12 ) =< S{fi)5{f2) > (2) 

Due to homogeueity and isotropy ^ depends only on the 
modulus of vector ri 2 . For the 3-point correlation function 
we have, in an analogous way, 

SP3 = p^\l+ ^(7-12) + C(^ 23 ) + ^(ris) + C,{ri2, r 23 , 7-13)] 

5Vi5V25V3 (3) 

where C(^i 2 T 23 ,J'ls) is the ’reduced’ 3-point correlation 
function or simply known as 3-point correlation function. 
The terms in eq. ^ depending on the 2-point correlation 
function reflect the excess number of triplets one has when 
there is an excess of pairs over a random distribution. Then 
the term C(>"i 2 , ^ 23 , ris) is simply the excess of triplets over a 
distribution with a given 2-point correlation function. When 
approximating by a continuous distribution of objects we 
have. 


C(7'i2,r23,ri3) =< S{fi)S(r2)S{r3) > ■ 


(4) 


Notice that we have just written < S{fi)S{f 2 )S{f 3 ) > instead 
of < 5 {fi)S{f 2 )S{f 3 ) >c, the cumulant or reduced moment, 
but because of the definition of 5{r} this one has a zero mean 
so both types of moments coincide. One interesting charac¬ 
teristic of the relation between moments and cumulants is 
the fact that for the case of a Gaussian distribution all the 
cumulants of order larger than 2 are zero [se e, for instance, 
(Stuart & Ord 1994; Bernardeau et al. 2001)]. 

In the study of the 3-point correlation function we will 
be interested in computing the Qs parameter in real space 


i n N-body simulations . The Q 3 parameter, introduced in 
( |Groth fc Peebles 1977 ), is defined as follows. 


Qs = 
Ch 


C{ri2,r23,ri3) 

CH{ri2,r23,ri3) 

= Gri2K{r23) + ?(n2)^(ri3) -k C(f’23)C(ri3). 


(5) 

( 6 ) 


where we have introduced a definition for the ’’hierarchical” 
3-point function (^h- The Q 3 parameter was thought to be 
a good descri ption of observations by being constant (see 
(Peebles 198C)), a result that is usually referred to as the 


hierarchical universe. As we will see below, in the weakly 
non-linear regime (or large scales compared to the correla¬ 
tion length), the universe is not exactly hierarchical as Q 3 
is not quite constant, but this is a good approximation in 
some cases. 


2.2 Tree-level contribution. 


In perturbation theory (starting from Gaussian initial con¬ 
ditions), the first non-null contribution to the 3-point func¬ 
tion is well-known and it can be even expressed in analytical 
form in the case of scale free (power law) power spectrum: 
P(k) = Ak^ [see (Peebles 198C; Fry 1984; Bernardeau et 
al. 2001 ) and references therein]. 


C(r-i2,r23,M) = 


10 n-l-3 ,ri 2 , r23, , 

— -I--!-) + 

17 n r23 ri 2 


(7) 


4(3-2(n-k3)-k(n-k3)V^) 

7 rP 


C(ri2)C(r23) + permutations 


C{ri 2 ,r 23 , E) = C(^i2T23, )"13(m)) where /r is the cosinus of 
the angle between ri 2 and r 23 , say a. An interesting prop¬ 
erty of the power law case is that C(^i 2 , >" 23 , ris) (and also 
Q 3 ) only depends on the ratios of distances ri 2 , r 23 , ris. As 
these distances bear the additional constraint of forming a 
triangle the dependence of g{ri2,r23, E reduced to only 2 
parameters, for example, u = r23/ri2 and v = (ri 3 —r 23 )/ri 2 
used by Peebles and collaborators, Jing and Borner and 
Buchalter and Kamionkowski or the parameters we will use 
here, u = r 23 /ri 2 and p = cos a where a is the angle de¬ 
fined above. We have chosen these latter parameters and not 
the former ones because they are in more direct relation to 
the way in which th e lis ts employed by our algorithm (ex¬ 
plained in subsection 3.2) are built. In figure we show how 
Q 3 {f'i 2 , r23,E) changes as a function of a for some values of 
u = r 23 /ri 2 = 1,2,10. The thicker lines correspond to a 
model with spectral index n = — 1 and the thinner ones to 
a model with n = — 2. It is noticeable how the Q 3 flattens 
with the decrease of the spectral index n. For each model 
parametrised by the spectral index there are three lines cor¬ 
responding to different u values. The solid ones are for u = 1 
(isosceles triangles), the short-dashed ones are for u = 2 and 
the long-dashed ones are for u — 10. 


2.2.1 CDM models 

In the case of CDM power spectra the solution is not ana¬ 
lytical but can be expressed in terms of some basic functions 
that can be easily computed numerically, 

Ciri2,r23,p) = y?(ri2)C(»'23) + 
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Figure 1. Plot of the value of Q 3 (ri 2 , ^ 23 ,/i) as a function a 
(recall that = cos a), the angle between ri 2 and ^ 23 , as obtained 
in perturbation theory (up to first non-null terms) for scale free 
models. Solid lines are for values of u = r 23 /ri 2 = 1, short- 
dashed ones are for u = 2 and long-dashed correspond to u = 
10. Different lines’ thicknesses correspond to different spectral 
indexes n as shown. 
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^{ri2)<j>'{r23) i{r23)(j)'{ri2) 


ri2r23 


r23 


ri2 


V(e(n2) + 3^^){ar23) + 3^^))- 

ri2 J'la 

-K^'{ri2)(l)'{r23) + {r23)<l>' iri 2 )) + permutations, (8) 

where 

.P(fc) sin(fcr) 


(j>{r) = I dk- 


/fc2 


kr 


(9) 


and P{k) = A k T(k)^ is the linear power spectrum for 
CDM models with primordial spectral index of the Harrison- 
Zeldovich type k, and ^ is a normalization factor. We take 
a transfer function T{k) of the Bond & Efstathiou type; 

T{k)=[l + {ak + {hkf-^+ {ckfY^ (10) 

with a = 6.4/rh“^Mpc, b = 3/rh“^Mpc, c = 

1.7/r/i“^Mpc, V = 1.13. r « Qrnh is the so-called ’shape 
parameter’ and h is Hubble’s constant in lGGKms~^ / Mpc 
units). This abo ve Eq.l^ is also identical to resu lts in (Jing 
& Borrer 1997; Gaztanaga & Bernardeau 1998) once a re¬ 


lation between their functions and ours is established. 

It is interesting to note that for CDM power spectra 
C,{r\2,r23, u) is no longer a function of the ratios of ri 2 , r23 
and ri 3 , as it was in the case of scale free power spectra, 
and so 3 variables are now needed (the extra variable de¬ 
fines an "effective” index in the analogy to the power-law 
case). In figure ^ we show this dependency on ri 2 by dis¬ 
playing two plots with the same range in the ordinate axis. 
On the top (bottom) panel we show the computed values 
for Q 3 in perturbation theory in the case ri 2 = 18/i“^Mpc 
(ri 2 = 12h“^Mpc). The solid lines correspond to isosceles 
triangles {u = 1) whereas long-dashed lines correspond to 
u = 2. The thick lines represent cosmological models with a 
shape parameter E = 0.2 and the thinner lines represent cos¬ 
mological models with F = 0.3. Notice how for small scales 


Figure 2. Plot of the value of Q 3 as a function of a, as obtained in 
perturbation theory (up to first non-null terms) for CDM models 
with shape parameter F = 0.2, thick lines, and F = 0.3, thin lines. 
Solid lines are for values of u = r 23 /ri 2 = 1 and long-dashed 
correspond to u = 2 . 


(the lower plot) the curves become flatter. It is also quite no¬ 
ticeable from figure ^ how increasing the shape parameter, 
F, mimics the effects of doubling the value of u for angles 
larger than, say, 70° but the difference is clear for smaller 
angles. In summary, both increasing the shape parameter 
and increasing scales make curves have more shape but the 
effects due to one case or the other can, in principle, be dis¬ 
entangled by looking at different ranges of angles above and 
below approximately 70°. 

Note the crossing of lines at a ~ 120° in both power-law 
and CDM cases, even for different values of u. The crossing 
occurs at a ~ 116 and Q 3 ~ 1.18, which seems to be a char¬ 
acteristic feature of the tree-level gravitational evolution. 

As we mentioned above this is the first non-null con¬ 
tribution (also called ’tree-level’) to the 3-point correlation 
function but there is an interesting feature in considering 
the next non-null terms, the so-called ’1-loop’ terms. 


2.2.2 Up to 1-loop contributions. 


Consider the second non-null terms in both C(r'i 2 , r 23 , p) and 
the denominator of Q3{r\2,r23, u) (the products of 2-point 
correlation functions), which were called Crr(ri 2 ,r 23 ,/r) in 
eq.^, the hierarchical 3-point correlation function, we have 
schematically, 


Q3 


(^(0) _|_ ^(1) 


( 11 ) 
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where is the first non-null term in perturbation the¬ 
ory of the denominator of Qs and the next non-null 
term. For the explicit conte nt of these terms in Fourier 
space see ( Scoccimarro 1997). For examp le, formulae from 
eq.[32] to eq.[36] in ( Scoccimarro 199^ give the expres¬ 
sion for the term in Fourier space (his B^^\ki,k 2 ,T)) 
and formulae eq.[39] and eq. [24-26] give the expression for 
the Fourier space counterpart of (his (fci, fc 2 , r)). 

and terms depend on the amplitude of the power 
spectrum squa red, A^, whereas and depend on 
(explicitly in (Scoccimarro 1997)). We will make use of this 
property in order to compare these terms extracted from 
simulations to the first non-null terms in perturbation the¬ 
ory iQf\rl2,r23,^^) = C™/Cl?^)- 


3 N-BODY RESULTS 

In this section we first describe briefly the N-body simu¬ 
lations that we will use. We next present two of the new 
contributions of this paper. The first of these contributions 
consists on an algorithm that drastically reduces the amount 
of time necessary to compute the correlation functions in 
N-bodies simulations with millions of particles. The sec¬ 
ond main idea refers to devising a way to disentangle non- 
linearities in the simulations by combining different output 
times. 


3.1 Description of the simulations 


The N-body simulations we will use here are evolved us¬ 
ing the P^M code of (Efstathiou et al. 1985). We have 
considered the following cases, A flat ACDM model with 
= 0.2 and F = 0.2 in a cube of side 378/i“^Mpc con¬ 
taining 126^ particles. We have also considered a SCDM 
(flm = 1, r = 0.5) model in a cube of side 378/i~^Mpc also 
containing 126® particles. We have 5 realisations for each 
model and consider two different output times as = 0.5 and 
0-8 = 1. Table |l| summarizes the main parameters in the 
simulations [see ^augh, Gaztanaga & Efstathiou 1995) for 


more details on these runs]. 

We have also use a new set of simulations with CDM 
shape F = 0.25 but with flm = 1 (called hCDM in table 
1^. These new simulations can be used to test the sensitiv¬ 
ity to Qrn with independent of the shape of P(k). They have 
a box size of 300h“^Mpc and twice the particle resolution 
of the previous sets {N = 200® particles). Thus they can 
also be used to check for shot-noise and resolution effects. 
The computation of gravitational dynamics was started at 
2 = 20 (instead of 2 = 10 in the other CDM sets), so that 
they are less sensitive to possible Zeldovich Approximation 
transients (see Baugh, Gaztanaga & Efstathiou 1995, Scocci¬ 
marro 1998). There are 5 realisations of the erg = 0.5 output. 


3.2 The algorithm to estimate correlations. 

Once simulations have evolved gravitationally from initial 
conditions we compute the 2-point and 3-point correlation 
functions. To do this with N particles, we would naively 
require and Y® operations for the 2-point and 3-point 
function respectively. The case of Y = 200® gives 5 • 10^° 



ACDM 

SCDM 

hSCDM 


0.2 

1 

1 

Qa 

0.8 

0 

0 

h 

0.7 

0.5 

0.3 

r 

0.2 

0.5 

0.25 

Npar 

126® 

126® 

200® 

Lsize 

378/i“^Mpc 

378/i“^Mpc 

300/i“^Mpc 


Table 1. Simulations’ main parameters. 


operations while the case of Y = 126® gives ~ 8 ■ lO*^® op¬ 
erations for the 3-point function. These numbers are out 
of reach for an ordinary computer if we bear in mind that 
many simulations (different realisations multiplied by out¬ 
put times) are n eeded. 

Based on (Frieman & Gaztanaga 1999) we devise a 
method which dramatically reduces the number of opera¬ 
tions. The first step is to discretize the simulation box into 
Lsize^ cubic cells. We assign each particle to a node of this 
new latticed box using the Nearest Grid Point particle as¬ 
signment (NGP). We thus obtain a density field over the 
Lsize^ lattice, resulting from just counting particles in each 
node. 

To compute now the 2-point and 3-point correlation 
functions in the lattice would naively take {Lsize^)^ and 
(Lsize^)^ operations respectively. This means a considerable 
reduction of operations when Lsize < Y^^® or, equivalently, 
when the mean cell density {N/Lsize^) is larger than 1. As 
we have seen in section both the 2-point and 3-point cor¬ 
relation functions depend only on lengths. This fact is used 
in the next subsection to devise a method which diminishes 
the number of operations even further. 


3.2.1 The 2-point correlation function. 

The 2-point correlation function depends only on the dis¬ 
tance between two points (labeled ’1’ and ’2’). In general, 
we want to consider a finite set of fixed distances which are 
much smaller than the box size (as otherwise our statistics 
will be dominated by sampling effects). So we would like to 
learn how to count only pairs of points separated by a fixed 
distance. It is then convenient to define a list of neighbors 
to each point. This list is formed by members of the lat¬ 
tice which are inside a spherical shell of a certain width of 
radius say ri 2 . In the lattice, this list is given by a single 
list of "displacement vectors” for all nodes. In a more gen¬ 
eral situation (where particles are not placed on a lattice) 
this list is different for each particle. In both cases the lists 
are pre-computed. In the case of the lattice, to calculate the 
list of "displacement vectors” takes a negligible amount of 
computing time. 

In order to count only these pairs, we run through all 
nodes in the cubic latticed box and for each node we loop 
over the attached list of neighbours. The saving in comput¬ 
ing time resides in the fact that for each node we do not 
compute the Lsize other nodes but only those nodes which 
belong to the list (say mi 2 ). The number of operations is 
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reduced from (Lsize^)^ to (Lsize^) ■ mi 2 , where mi 2 (the 
number of the nodes belonging to a spherical shell) is always 
much smaller than Lsize^ (the nodes in the box; the diam¬ 
eter of the shell is restricted to be smaller than the side of 
the box). To show some numbers, let us consider the case 
with Lsize = 63. Naively the number of operations would 
be (63®)^ « 6.2 • 10^°. When we are interested in a hxed 
scale, say, for example ri 2 = 3 (in lattice units) there are 66 
entries (i.e. nodes) in the list of ’’displacement vectors”, so 
the number of operations would be (63®) • 66 « 1.6 • 10^. An 
additional saving of a factor of two can be achieved by realis¬ 
ing that the ordering of the pairs is not relevant. This means 
in practice that we only need a semi-sphere of neighbours in 
our list. 

One might think that this saving is ficticious as it would 
take the same time to loop over all particle pairs as to loop 
over all possible pair distances using the neighbours’ list. 
But recall that in practice the scales of interest are always 
ri 2 « Lsize, so that it does not make sense to loop over 
all pair distances 

The advantage of this method in computing time be¬ 
comes far more evident in the case for the 3-point and higher 
order correlation functions. 

3.2.8 The 3-point correlation function. 

The 3-point function depends only on the length of the three 
vectors formed by the 3 point-particles, that is on the sides 
of a triangle formed by these 3 point-particles being at the 
vertexes. We can also express this dependency as a depen¬ 
dency on any two vectors, chosen from the three available, 
and the angle between them a. Then the method now con¬ 
sists of using the same lists of neighbours, as for the 2-point 
function’s case, with two nested spheres, each one of radius 
equal to the length of one of these two vectors. Thus, we run 
through every single lattice node in our simulations boxes 
(node labeled ’1’ in figure |^. To each node in the lattice 
we run over the list of neighbours within a semi-sphere shell 
of radius ri 2 and a certain width (all these nodes are candi¬ 
dates for being labeled as ’2’ in figure |^. Each of the new ’2’ 
nodes of the semi-sphere are now centers for another spher¬ 
ical shell of radius r 23 and a certain width (these new nodes 
are labeled ’3’ in figure In summary, the list of triplets 
separated by distances ri 2 and r 23 is built from two nested 
list of neighbours of separations ri 2 and r 23 (see figure |^. 
As in the case of the 2-point function, the first semi-sphere 
(as opposed to a full sphere) is used to avoid double counting 
pairs. 

This method of building the neighbours list associated 
to a node in the box makes us save even more time than in 
the case of the 2-point correlation function. The number of 
operations to be done in the case for fixed scales are reduced 
to (Lsize^) ■ mi2/2 • 17123 instead of (Lsize^)^, where mi2/2 • 
17123 is the number of nodes in the neighbours list. To write 
down some numbers let us take the example of latticing the 
box with Lsize = 63. To compute the 3-point correlation 
function naively would take (63®)® ~ 1.6 • 10^® operations 
in the lattice (recall that this was as large as 5 • 10®® from 
the raw particle positions) whereas by using our method of 
associating a neighbours list for hxed scales, say for instance 
ri 2 = 723 = 3 (in lattice units), would take (63®) • 33 • 66 « 
5.4 • 10® operations. 



Figure 3. Schematic plot for the method explained in the text 
to compute the 3-point correlation function in simulations. From 
point labeled ’1’ we trace the shell of a semi-sphere of a finite 
thickness with radius ri 2 . From each point of the shell (points 
labeled ’2’) we trace another spherical shell with radius r 23 (a 
whole sphere this time). 


This method can be trivially extended to higher order 
correlations. 


3.2.3 Generalities of the algorithm. 

One characteristic of our algorithm is that shells, from which 
we build our lists, are distorted. This is always the case 
when trying to obtain spherical objects in a cubicly grid- 
ded box. Obviously, the distortion in the spherical nature of 
shells is going to diminish with the increase of the radius in 
terms of how many cells long the latter is. We have stud¬ 
ied these effects by using different size cells (or lattices) to 
study the same physical side. For example, in the case of 
ACDM model with Qm = 0.2 simulation, which has a box 
of 378/i~^Mpc, we have used a lattice of Lsize = 63 nodes 
with a radii ri 2 = 3 and r 23 = 3 in cells units (to study the 
physical scales ri 2 = r 23 = 18/i“^Mpc). We also latticed the 
simulation with with Lsize = 126 and a radii ri 2 = r 23 = 6 
(in cell units), which corresponds to the same physical dis¬ 
tance. We find that the differences are negligible when we 
use 3 or more cell units for the distance. To reduce the num¬ 
ber of operations we want Lsize to be as small as possible, 
we therefore choose Lsize so that the scale of interest corre¬ 
sponds to r = 3 cell units. The thickness of the shell is also 
a parameter which affects the distortion, as some cubic cells 
may be missed out if we choose this parameter to be too 
small and too many cells will be included if this parameter 
is too large. We have studied these effects and have adopted 
a compromise of a thickness of 20% in cell side units. 
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3.3 Disentangling nonlinearities in the 2-point 
and 3-point correlation function. 


As can be seen from its definition Q 3 depends on two sides 
of the triangle formed by the three points and on the angle 
between them. Examining the a dependence it is clear that 
nonlinearities are going to be important when a ^ 0. When 
the third side approaches zero the 2 -point correlation func¬ 
tion on this scale is going highly non-linear. On the other ex¬ 
treme, at large scales, it is also possible that Q 3 diverges be¬ 
cause the denominator approaches zero. This is unavoidable 
in CDM models, as the 2-point correlation becomes negative 
at some large scale which is not very far from the non-linear 
scale. This divergence of Q 3 has promoted proposals for ’bet¬ 
ter behaved’ normalizatio ns for C. but with similar scaling 
as Q 3 . This is the case of (Buchalter & Kamionkowski 199£) 
where they define = C(’'i2,^’ 23 ,m)/o' 8 with: 


= 6 (ri 2 = 0 ,Rf= 8 ) = J dkP{k)W{kRf ( 12 ) 

and W (kR) is the Fourier transform of the window function. 
Non-linearities affects both the numerator and denominator 
of Q 3 . We need to separate the different non-linear compo¬ 
nents in Q 3 if we want a consistent comparison of PT with 
simulations. In order to estimate the effects of mild nonlin¬ 
earities in the simulations we devise the following method, 
based on the scaling of both the 3-point and 2-point cor¬ 
relation functions with the amplitude of the linear power 
spectrum (or growth factor) in perturbation theory. 

As explained in the previous section the tree-level (that 
is, the first non-zero term) contribution to the 3-point func¬ 
tion depends on the square of the amplitude of the power 
spectrum and the one-loop level (next order in perturbation 
theory) depends on A?. We can do exactly the same in the 
denominator developing it up to order A®. Imagine now that 
we get two outputs of our simulations with different ampli¬ 
tudes, i.e. different erg. For example, one with ug = 0.5 and 
another one with erg = 1. This is a difference in amplitude 
of a factor of 4. We can then pose a very simple system of 
two equations with two unknowns, 


X + y = ni 
X y 

16 M 


n2 


(13) 


where x stands for or (i.e., tree-level terms), y stands 
for or (i.e., 1 -loop terms) and ni and n 2 are the 
results we obtain for the 3-point function with simulation 
output CTg = 1 and ug = 0.5 respectively. So, in this way we 
are able to obtain the tree-level (’x’) both in the numerator 
and denominator. 

This could obviously be extended to higher order cor¬ 
rections if we had more output times. Note that these set 
of equations will always return a solution. The solution will 
only be meaningful if deviations from the tree-level are small 
enough at the order of perturbations we are considering. One 
way to test the consistency of the solution is to verify that we 
are able to recover the tree-level solution (’x’) from the dif¬ 
ferent outputs. Clearly, on small, strongly non-linear scales, 
the perturbative approach will break down and the solution 
will be meaningless (in such case we do not expect to recover 
the tree-level solution for ’x’). 



a 


Figure 4. Plot of the value of Q 3 (ri 2 , ^ 23 , o) as a function a 
in ACDM models with erg = 0.5, the angle between ri 2 and 
^ 23 , for different values of |ri 2 | = 1^231 = 24,18,12,6/i“^Mpc 
(opened circles, closed squares, open squares and triangles). The 
corresponding lines (continuous, dotted, short-dashed and long- 
dashed) show the tree-level perturbation theory. 


4 COMPARISON TO SIMULATIONS 


In this section we compare the results to actual simulation 
values. We first analyze different triangles configurations 
(isosceles and equ ilateral) using our algorithm (explained 
in subsection |3.2|) and examine how they compare to theo¬ 
retical values in subsections 0 and [ 1.4 In subsection |4.3| 
we apply our method to disentangle possible ’non-linearities’ 
due to 1 -loop terms and ob tain a the tree-level from simu¬ 
lations (as described in 3.3). 


4.1 Shape dependence for isosceles 

Figures ^ and ^ show the estimations of Q 3 in the 
ACDM model for isosceles triangles |fi 2 | = IfisI = 
24, 18,12, 6 /i“^Mpc at two different output times: erg = 0.5 
in fig. and (Jg = 1.0 in fig. The figures with error-bars 
show the mean and variance in 5 simulations, while the dif¬ 
ferent lines show the corresponding leading order (tree-level) 
PT prediction. As can be seen in fig|^ the tree-level is re¬ 
covered well for all scales |fi 2 | = |r 23 | > 6 /i~'^Mpc/h in the 
erg = 0.5 output. At |fi 2 | = |r 23 | = 6 /i~^Mpc//i there are 
significant deviations from the leading order results indicat¬ 
ing that loop corrections are important. 

In the erg = 1.0 output (fig.^) the discrepancies are 
significant even for jfi 2 | = |r 23 | = 24h“'^Mpc. Note how 
deviations from tree-level tend to smooth or wash out the 
shape dependence for |fi 2 | = If^sj > 12/i“^Mpc, while on 
smaller scales there is more shape dependence than in the 
tree-level, in the sense that Q 3 is smaller at smaller a and 
larger at larger a. This tendency is also true when com- 
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Figure 5. Same as fig.k 


for the (Tg = 1 output. 


paring the shape of two outputs at a given scale: i.e. the 
later outputs have more shape dependence than the earlier 
outputs at small scales (< 12/i~^Mpc), while the opposite 
tendency occurs for the larger scales. This effect looks sim¬ 
ilar to th e " previrialization effect” pre s ent in the 2 -poin t 
function ( Scoccimarro fc Frieman 1996| ; Lokas et al. 1996 h 
A simil ar tr end was also found in the bispectrum (Scocci- 
marro ^ 997). The idea is that there is a change of sign in the 
1 -loop contribution which, for a power-law initial spectrum, 
is related to a critical value in the spectral index. We will 
study this in more detail below. 


Figure 6. The normalized 3-point funciton gg for equilateral tri¬ 
angles as a funciton of the triangle side. The continuous line cor¬ 
respond to the PT tree-level predictions while the triangles and 
errors show the corresponding estimation for the mean and vari¬ 
ance of 5 realisations in the eg = 1.0 output. Opened squares show 
the mean of the 5 realisations in the erg = 0.5 output. The top 
panel correspond to the SCDM model while the bottom panel 
shows the ACDM results. 


—2.2) PT corrections on Qg are below tree-level. On smaller 
scales, i.e. ri 2 < 4/i“^Mpc, non-linear effects take over and 
Qs increases again above the PT results. These deviations 
from the leading order results are small, but significant given 
the sampling errors. 


4.2 Equilateral configurations 

Fig.p shows Qg for equilateral triangles (i.e. fixed |ri 2 | = 
|r 2 g| and a = 60°) as a function of the triangle side 
(|ri 2 | = l^ 2 g|), for both the erg = 1.0 output (triangles) 
and the erg = 0.5 (squares), in both the SCDM (top panel) 
and the ACDM (bottom panel) model. The dashed line just 
connects the points at different scales. Also shown as a con¬ 
tinuous line is the leading order (tree-level) contribution. 

For both models, the erg = 0.5 output (squares) shows 
good agreement with the PT results. For erg = 1.0 the 
SCDM shows significant deviations from PT at all scales. 
On the largest scales the 2-point function goes to zero at 
smaller scales in the SCDM model which explains why the 
errors are much larger around 20/i“^Mpc in this model. 

The ACDM model shows larger Qg values than PT at 
ri 2 — 20/i“^Mpc and smaller values at ri 2 — 6 /i“^Mpc. This 
is related to the change in shape away from the leading order 
result noted in the subsection above (e.g. see figj^. More in 
detail, at large weakly non-linear scales, ri 2 > 12h“^Mpc 
where the 2 -point function is steeper (with a slope 7 < — 2 . 2 ) 
PT corrections on Qg are above tree-level, while for ri 2 < 
12h“^Mpc the 2-point function is flatter (with a slope 7 > 


4.3 Loop corrections 

In figure we show our method to recover the tree-level and 
the first-loop contributions for the denominator and numer¬ 
ator of Qg, i.e . for (■, and C,h respectively, as explained in 
subsection 3.3. All the plots in figure ^ are for a flat ACDM 
model with Q.m =0.2. The plots on the left column are for 
isosceles triangles with equal side of 18/i“^Mpc long whereas 
the ones on the right correspond to isosceles triangles with 
equal side of 9/i“^Mpc long. From top to bottom, we show g 
(the 3-point function), C^h (the hierarchical combination in 
eq. Q) and Qg (the ratio of both). In all the plots of figure 
Q the solid lines correspond to the tree-level theoretical pre¬ 
dictions from PT, the long-dashed and short-dashed lines 
correspond to the tree-level and 1 -loop (ng = 1 ) estimation 
from the 2 simulation output times, as described in eq. |^. 

We can now address the question raised in the previous 
subsections of whether or not the relative change of sign in 
the deviations of Qg from tree-level PT (e.g. shown in fig.^ 
and figj^ is related to the ’’previlarization effect” on the 
2 -point function. 

In the middle row plots we can see how for both scales 
(18/i“^Mpc and 9/i“^Mpc) the 1-loop contribution to (^h 
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(short-dashed line) is lower than the tree-level one (long- 
dashed line); this can be interpreted in both cases as the 
well-know n effect of previrialization found in the 2-poin t 
function ( gcoccimarro fc Frieman 1996| ; Lokas et al. 1996). 
But here, we have to bear in mind that (^h is a product of 
two 2-point correlation functions. 

In the top panels one can see how the 1-loop contribu¬ 
tion to (short-dashed lines), are larger than the tree-level 
contribution around a ~ 60 deg. for 18/i“^Mpc, while they 
are smaller than tree-level for 9h“^Mpc. This effect is the 
one mostly responsible for the relative change shown in Qs 
(smoothing/enhancing the shape with respect to the tree- 
level). Thus, we see that this tendency is not just given by 
the same critical scale index in the 2-point previrialization, 
as here we found a case where there is change of sign at 
the 3-point function level (when comparing 18/i“^Mpc with 
9h“^Mpc results), but not for the 2-point function. This is 
not surprising, as there are different PT terms contributing 
to each of these expressions. 

In figure we compare the 3-point correlation func¬ 
tion {(), the hierarchical 3-point correlation function (Crr) 
and the full Qs parameter for the model hSCDM, on the 
left column, which has Qm — 1, h = 0.3 (see table to 
the corresponding values in the SC DM model, on the right 
column, with = 1, h = 0.5 (the old standard, see table 
^ ). As before the solid lines correspond to the theoretical 
results in PT. The short-dashed shows the us = 0.5 out¬ 
put. On the left panel, the short-dashed and long-dashed 
lines correspond to the tree-level and 1-loop reconstruction 
in the SCDM case. The purpose of this figure is to illus¬ 
trate how the C,H and Qa functions change with the cos¬ 
mological model. The hSCDM model has an approximate 
shape parameter for the power spectrum P = 0.25 whereas 
the SCDM model has a P = 0.5. Notice that although the 
shape of the curves of the left and right panels in hgure ^ 
are quite similar the actual values in the range of the plots 
are different. 

Figure ^ also illustrates how 1-loop corrections (shown 
as short-dashed lines) in the SCDM tend to flatten the Qa 
curve (see bottom right plot), in a similar, but more dra¬ 
matic, way as we found in the ACDM model in hgure 


4.4 Other configurations 

In hgure ^ we show plots for Qa for the ACDM model with 
u = r2a/ri2 = 1 (top row), u = 2 (middle row) and u = 3 
(bottom row) for two values of ri2: ri2 = 18h“^Mpc for 
the left column and ri2 = 12h“^Mpc on the right one. The 
solid lines correspond to theoretical perturbation theory, the 
long-dashed lines correspond to the tree-level recovered from 
simulation by using our method in eq. the short-dashed 
lines correspond to the 1-loop level {as = 1) recovered from 
simulations. Note the diherence range of values of Qs in each 
panel. The purpose of this hgure is to compare the effects 
at different scale lengths. On the top left plot one can see 
how at somewhat large scales, a > 100 deg. simulations 
are below the theoretical prediction. This ehect is due to 
volume effects (or small statistical sampling); notice that 
a > 100 deg corresponds to ris > 43h~^Mpc which is above 
one tenth of our box side, 378h“^Mpc. This type of volume 
effects produce a (ratio) bias in the estimators which is not 
taken into account by averaging the results over different 


realisations ( Hui fc Gaztanaga 1998| ). For larger u this effect 
is even worse. This can be seen in the bottom left panel, 
where the error bars increase and the estimation becomes 
quite uncertain. 

At smaller scales (right panels) volume ehects are 
smaller and there is a general agreement of tree-level PT 
with simulations. Higher order corrections are small. 


5 CONCLUSIONS 


In this paper we have studied how to obtain the 3-point cor¬ 
relation function from simulations and we have compared 
them to the predictions of perturbation theory (PT). In 
order to obtain the 3-point correlation function from sim¬ 
ulations in a reasonable amount of time we have devised 
a method which dramatically reduces the number of op¬ 
erations to be carr ied out. The method, explained in de¬ 
tail in section 3.2, basically consists of assigning a list of 
neighbours to every node of a previously gridded box. This 
method, allows us to compute the i^, C,h and Qs for several 
million particle s in only a few minutes. 

In section 3.3 we also propose a new method to sepa¬ 


rate from simulations the tree-level and 1-loop contribution 
to the 2 and 3-point functions. This allows a more clear com¬ 
parison to the theoretical predictions. We obtain that on suf¬ 
ficiently large scales, simulations results are in good agree¬ 
ment with leading order perturbat ion theory. This seems in 
contradiction with the results of (Matsubara & Suto 1994; 
ling & Borner 1997); but the contradiction is only apparent. 


In these previous analysis they focus on the results of small 
scales (less than 8/i~^Mpc). Here we have shown that next 
to leading order (i.e. 1-loop) corrections are important on 
these scales and that one needs to consider larger scales (or 
earlier times) in order to reach the leading order (tree-level) 
regime. Once this is taken into account there does not seem 
to be a contradiction with wha t was found in (Matsubara & 
Suto 1994; Jing & Borner 1997). On even larger scales, the 2- 
point correlation goes to zero and the results are dominated 
by sampling errors. The goodness of the agreement and the 
specihc scales where PT is valid, depend on the shape of the 
initial fluctuations (i.e. the P parameter in CDM models), 
the redshift (as), the particular 3-point configuration under 
study and the level of accuracy in the comparison (i.e. the 
sampling errors). For the currently favoured ACDM model 
at (78 = 1 there is reasonable agreement at leading order for 
scales larger than 20/i“^Mpc. 

We also characterise next to leading order (1-loop) cor¬ 
rections and find that in the ACDM model they seem 
to change sign around 12/i“^Mpc in Qs (see fig.^). At 
12h~^Mpc, loop corrections seem to cancel and the tree- 
level PT prediction gives a better match to the simulations 
at this scale than at slightly larger scales (this statement 
would depend on the accuracy of our determination, e.g. 
the si ze o f the sample, as errors increase with scale). In sec¬ 
tion §4.3 we hnd that this change of sign in the 1-loop term 


is not clearly correlated to the cri tical index found in the 
so ca l led ’’previrializatio n effect” (|Scoccimarro fc Frieman 


1991 ; Lokas et al. 1996 ) f or the 2-point fun ction and also 


found in the bispectrum (Scoccimarro 1997). 

We have noted a crossing of the different shape predic¬ 
tions at Q ~ 120 degrees in both power-law and CDM cases. 
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0 50 100 150 0 50 100 150 



0 50 100 150 0 50 100 150 


a a 

Figure 7. Recovering the tree-level and 1-loop contributions to the numerator (top) and denominator (middle) of Qz, and also for Qz 
itself (bottom). The cosmological model is flat ACDM with Qrn = 0.2. The left (right) column correspond to isosceles triangles with 
equal side being 18/i~^Mpc (9/i“^Mpc) long. The solid lines represent tree-level PT prediction, the long-dashed lines shows tree-level 
recovered for the simulations while the short-dashed corresponds to the estimated 1-loop level in the simulations (at as = 1). For clarity, 
only the larger error-bars in either of the two dashed lines are shown. 
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r,2=18Mpc/h, u = 1, hSCDM r,2=18Mpc/h, u=1, SCDM 




r,2=18Mpc/h, u = 1, hSCDM 



r,2=18Mpc/h, u=1, SCDM 



r,2=18Mpc/h, u = 1, hSCDM ri2=18Mpc/h, u=1, SCDM 




a 
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Figure 8. On the left column we show results for s,nd Qz in the hSCDM model with Qrn = 1 h = 0.3, corresponding to isosceles 

triangles configurations with equal side 18h“^Mpc. The dashed line with error-bars correspond to the as = 0.5 output, the continuous 
line shows the leading order PT prediction. On the rmht column we show the same case for the model SCDM (Qm = l,h = 0.5 ). Here 
the criteria for lines are identical to those for figure M. 
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r,2=18Mpc/h, u = 2, ACDM 
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ri2=18Mpc/h, u = 3, ACDM 
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ri 2 = 12 Mpc/h, u = 2, ACDM 
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ri 2 = 12 Mpc/h, u = 3, ACDM 
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Figure 9. On the left column we show results for the ACDM model with ri 2 = 18h~^Mpc and on the right we show results for the 
same model with ri 2 = 12h~^Mpc .The solid lines correspond to theoretical perturbation theory, the long-dashed lines correspond to 
the tree-level recovered from simulation by using our method in eq. O and the short-dashed lines correspond to the up to 1-loop level 
(78 = 1 recovered from simulations. We only show error bars in the curve in which they are larger. 
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even for different values of u. That, is when one of the in¬ 
terior angle of the triplet is about 120 the value of Q 3 is 
Q 3 ~ 1.2, independently of the scales involved!. This seems 
to be a characteristic feature of the tree-level gravitational 
evolution, but only holds approximately after loop correc¬ 
tions. We do not yet have a physical interpretation for this 
feature, but it could be used to test gravitational instability 
against observational data. 

Surpris ely, even though PT ha s been known for many 
years now (Peebles 198C; Fry 1984), the results presented 


here are the hrst ones to show a comparison of PT with 
the 3-point function in N-body simulations at the appro¬ 
priate scales. Most of previous comparisons of the 3-point 
s tatistics in PT to simulations were done in Fouri e r space 
(Scoccimarro fc Frieman 1996; Bcoccimarro 1997; Scocci- 


marro jOPCj ) , in angula r space ( FViemau&_GnztahngaM99f ) 

or for smoothed fields (BernaxdeauM994 ; [Baugh, Gaztanaga 
& Efsti thiou 199^ ; [Bernardeau et al. 200l| ). The compar- 
ison presented here for the 3-point function (in configura¬ 
tion space) extend, for the first time, to the weakly non¬ 
linear s cales were this compari son is meaningful. Not Sur¬ 
prisely (Bernardeau et al. 2001), we find an excellent agree¬ 
ment with tree-level PT theory. Moreover we characterise 
the range of validity of this agreement and the next order 
(1-loop) corrections. 
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